2024.3.4 ４次のルンゲクッタ法

$ x_{n+1} = x_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)

$ t_{n+1} = t_n + h

ただし、

$ k_1 = f(t_n, x_n)

$ k_2 = f(t_n + \frac{h}{2}, x_n + \frac{h}{2}k_1)

$ k_3 = f(t_n + \frac{h}{2}, x_n + \frac{h}{2}k_2)

$ k_4 = f(t_n + h, x_n + h k_3)

code:rk4.py

import numpy as np

import matplotlib.pyplot as plt

def func_xd(t, x):

return A@x

def solve_rk4(times, x, h):

x_hist = []

for t in times:

x_hist.append(x.flatten())

k1 = func_xd(t, x)

k2 = func_xd(t + h/2, x + k1*h/2)

k3 = func_xd(t + h/2, x + k2*h/2)

k4 = func_xd(t + h, x + k3*h)

x = x + h*(k1 + 2*k2 + 2*k3 + k4)/6

return np.stack(x_hist)

##

t_initial = 0 # 初期時刻

t_final = 10 # 終端時刻

t_interval = 0.01 # 時間の刻み幅

m, c, k = 1.0, 1.0, 1.0 # 質量、減衰係数、ばね定数

##

x_initial = np.array(x_initial)

times = np.arange(t_initial, t_final, t_interval)

result = solve_rk4(times, x_initial, t_interval)

plt.grid()

plt.legend()

plt.show()